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ABSTRACT 

We present a new method to identify large scale filaments and apply it to a cosmolog- 
ical simulation. Using positions of haloes above a given mass as node tracers, we look 
for filaments between them using the positions and masses of all the remaining dark- 
matter haloes. In order to detect a filament, the first step consists in the construction 
of a backbone linking two nodes, which is given by a skeleton-like path connecting 
the highest local dark matter (DM) density traced by non-node haloes. The filament 
quality is defined by a density and gap parameters characterising its skeleton, and 
filament members are selected by their binding energy in the plane perpendicular to 
the filament. This membership condition is associated to characteristic orbital times; 
however if one assumes a fixed orbital timescale for all the filaments, the resulting 
filament properties show only marginal changes, indicating that the use of dynami- 
cal information is not critical for the method. We test the method in the simulation 
using massive haloes(M > lO 14 h _1 A/ ) as filament nodes. The main properties of 
the resulting high-quality filaments (which corresponds to ~ 33% of the detected fil- 
aments) are, i) their lengths cover a wide range of values of up to 150h _1 Mpc, but 
are mostly concentrated below 50h~ 1 Mpc; ii) their distribution of thickness peaks at 
d = 3.0h _1 Mpc and increases slightly with the filament length; iii) their nodes are 
connected on average to 1.87 ± 0.18 filaments for ~ 10 141 Af Q nodes; this number 
increases with the node mass to ~ 2.49 ± 0.28 filaments for ~ 10 14 9 M© nodes; iv) on 
average, the central density along the filaments starts at almost a hundred times the 
average density in the regions surrounding the nodes and then drops to about a few 
times the mean density at larger distances, where it remains roughly constant over 
20 to 80% of the filament length (this result may depend on the filament length); v) 
there is a strong relation between length, quality and how straight a filament is, where 
shorter filaments are those characterised by higher qualities and more straight-line like 
geometries. 

Key words: large scale structure of Universe. 



1 INTRODUCTION 

The large scale distribution of galaxies and dark mat- 
ter(DM) shows a web-like structure composed by clusters, 
walls, filaments and void regions, and is usually referred to as 
the cosmic web. These structures can be easily detected by 
eye in numerical DM simulations or in the observed distri- 
bution of ga laxies in large surv eys such as the Sloan Digital 
Sky Survey (|York et al. Il2000l . SDSS). 

For clusters and voids, there are several well estab- 
lished automated identification methods which have been 
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broadly used, such as the Friend-of-Frien ds algorithm for 
halo/cluster detection ([Davis et al. I Il985l . FOF), and the 
IPadilla. Ceccarelli fc Lambas "(12005 ) algorit hm for reliable 
detection of voids (see Icolberg et al. I [2008 for a complete 
review on different void detection methods). In the case of 
filaments and walls this task is markedly difficult since, in 
general, there is still no clear consensus on how to charac- 
terise them; filaments and walls show complex 3D shapes. 



There are different approaches to the study of filaments. 
From the theoretical point of view it was found that the 
gravitational collapse of matter on lar ge scales leads t o 
the formation of s heets and filaments (jZel'dovich I Il970l ) . 
iBond et al. I §996) studied tidal fields in the large-scale 
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structure (LSS) and showed how these produce filamentary 
structures. 

There are several sets of filaments which have been 
identified an d characterised by eye in both simulati ons and 
observations. IColberg. Krughoff fc Connolly I (|2005l ) identi- 
fied by eye 228 filaments between massive neighbouring 
haloes in a DM simulation, and described several interest- 
ing statistical properti es us ing this sample. In observations, 
iPimbblet et al. I (|2004h and lPorter et all (120081 ) identified fil- 
aments in large surveys by eye, and dark matter (DM) fil- 
aments were also be detected between clusters o f galaxies 
using weak lensing techniques l|Mead et al. I 2002). In x-ray 
observations, it has also bee n possible to detec t hot gas fil- 
aments connecting clusters (|Scharf et al. |[2000l ) . 

The study of statistics and the topology of the galaxy 
distribution with the aim to search for filaments starts 



very e a rly, with studies by Zel'd ovich. Einasto & Shandarin 
19821 ) JShandarin fc Zel'dovich I (|l983l )T and lEinasto et al. I 
19841 ). Options to automate the search of filaments in- 



clude the use of statistics on the morphology of struc- 
tures, such as Minkowski functionals, minimal spanning 
trees (M ST), percolation method s and shapefinders (see re- 
view by iMartmez fc Saarl I2002T I . The mini mum spanning 
tree m ethod was introduced in cosmology by I Barrow et al.l 
(1985). This produces a unique graph which connects points 
of a process without closed loops, but describes mainly the 
local nearest-neighbour distribution and is una ble to provide 
a full characterisation of the LSS. Shapefinders (jSahni et al. I 
1 199ST > have also been used to identify filaments. 

In three dimensions, the morphology of a compact man- 
ifold can be characterised by four Minkowski functionals: 
volume, surface area, integrated mean curvature and inte- 
grated gaussian curvature. It is possible to define a number 
of quantities related to those functionals; if a set of positions 
of galaxies or haloes is characterised by particular values of 
ratios between the Minkowski functional s, it is very likely 
that it will show a filamentary shape (|f3haradwai et alTl 
2000), but this does not guarantee a true detection of a 
filament or that all the selected members actually belong to 
the filament. 

Anothe r algorithm for t he detection of filaments was 
proposed bv lPimbbletl (|2005l ) based on the assumption that 
the orientations of constituent galaxies along such filaments 
are non-isotropic. This method works well on straight fila- 
ments with separations smaller than 15Mpc/h, as has been 
shown in their application to the 2-degree Field Galaxy Red- 
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shift Survey (2dFGRS. IColless et al 1120011). 



The Skeleton method l|Eriksen et al. I 12004 
iNovikov etaTI 120061 ) has proven useful for the detec- 
tion of possible filamentary structures in continuous two 
dimensional density fields. The skeleton is determined by 
segments parallel to the gradient of the field connecting 
saddle points to local maxima. The method involves 
interpolation and smoothing of the point distribution, 
introducing the kernel band- width as an extra parameter in 
the procedure of estimatin g the density fi e ld. Ex tending this 
work to three dimensions, Sousbie et al. I (|2008l ) found good 
agreement between detected skeletons and eye detections 
in a numerical DM simulation. By using the Hessian 
matrix eigenvalues t hey were able to de t ect fil amentary 
structures (See also Aragon-Calvo et al. 1 l2007al . 2007b). 
iBond. Strauss fc Cen I (|2009l) also use the Hessian matrix 



of the galaxy density field smoothed on different scales to 
characterise the m orphology of the LSS in mock catalogues 
and in the SDSS (|Stoughton et al. 1120021 ); they use their 
detected structures to determine the typical scales where 
filaments, clumps and walls are dominant. 

The Candy Model used by IStoica et al. I (|2005l ). is a 
two-dimensional marked point process where segments serve 
as marks. This method has been adapted to three dimen- 
sions and also impr oved to a more general Bisous Model 
|Stoica et al. I [2008), producing detections in very good 
agreement with the result of eye detection in tracing fila- 
mentary structures using only galaxy positions (as in the 
method we will present). However, the detection and thick- 
ness of the resulting filaments is only given by a coverage 
threshold (percent of total points, to be included in fila- 
ments). 

The sp in and orientation of haloes i n filamen ts has been 
studie d by lAragon-Calvo et al. I (|2007bl ) and IZhang et al. I 
(2009). They use a Multi-scale Morphology Filter (MMF) 
and compute the Hessian Matrix eigenvalues in a density 
field smoothed on different scales, to divide the full volume 
of their samples into cluster, filament and wall like struc- 
tures. However, this method, as well as other Hessian matrix 
based methods, is affected by a lack of an ability to deter- 
mine the thickness of filaments, and are difficult to apply 
to observational data, where one needs to define whether a 
galaxy is a member of a cluster, filament or void. 

In this work we propose a new automated method 
to detect filaments which builds upon ideas of several of 
the methods mentioned previously. A novel feature of the 
method is that it is designed to search for filaments us- 
ing nodes (corresponding to haloes o r galaxy clusters as in 
IColberg. Krughoff fc Connolly I l2005l ) selected by applying 
lower limits on their mass (or proxy for mass). This new 
method aims to be applicable to discrete halo or galaxy po- 
sitions even when these are so sparsely distributed that it is 
not possible to define a smooth density field, or that the Hes- 
sian matrix cannot be computed with an adequately high ac- 
curacy. This makes it particularly suitable for observational 
data such as the 2dFGRS or SDSS. In addition, we replace 
the smoothing scales and filament coverage thresholds by 
parameters with improved physical meaning. In this new ap- 
proach a filament quality depends on parameters related to 
the relative density and gaps of the filament skeleton, and its 
members are identified as the haloes or galaxies with binding 
energies with respect to the filament in the plane perpendic- 
ular to its skeleton. We will use the numerical simulations to 
calibrate the binding condition using objects with a collapse 
time and radius that can be computed even when dynami- 
cal information is not available, as is usually the case with 
observational data. In the latter, measurements or proxies 
for galaxy masses will still be required in order to define the 
filament membership condition. 

This paper is organised as follows. Section 2 presents the 
numerical simulation on which we perform our automated 
search for filaments. The method is presented in Section 
3, which also includes details on the measurement of the 
local density field, and describes the input parameters of 
the algorithm. Section 4 presents the results and Section 5 
concludes this work with our conclusions. 
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2 THE NUMERICAL SIMULATION 

We use a cosmological DM simulation with parameters cor- 
responding to the concordance ACDM model (cold dark- 
matter, n b = 0.045, n DM = 0.235, Q DE = 0.72, h = 
0.72, a s = 0.847, & n = 1), 500 3 particles and a peri- 
odic cube side of 250Mpc/h. At z = we find 176,041 
haloes and subhaloes in the mass range 1.4 x 10 11 h- 1 M Q < 
M < 1.5 x lO^h" 1 ^, i dentified using the AHF code 
l|Knollmann fc Knebe ll2009t) . For the detection of filaments, 
we select as nodes a total of 427 haloes with M > 
10 14 1i -1 Mq. The node pairs that will be the candidates 
for filament search are constructed using neighbour nodes, 
which are easily obtained using Voronoi Tessellations(VT 
hereafter, to be explained in more detail in the next sec- 
tion). We obtain a total of 3,385 node pairs with separa- 
tions < 65h~ 1 Mpc, using periodic conditions (310 node pairs 
straddle the simulation borders) ; Figure[T]shows all the node 
pairs in a slice of the simulation. In the next section we will 
apply the filament detection method to each of these node 
pairs. 



3 METHOD 

Our filament detection method is described in this section. 
We apply the method to dark-matter halo positions in the 
simulation as a first step towards the detection of galaxy 
filaments from observational data. A future extension will 
also use halo substructure as well as galaxies from a semi- 
analytic model so as to mimick real galaxies as closely as 
possible (as galaxies are thought to form in the potential well 
of DM haloes and subhaloes). When applying our method to 
semi-analytic galaxies we will be able to detect the effects of 
using proxies for the host dark-matter halo masses obtained 
from a galaxy catalogue (e.g. dynamical masses, luminosities 
in different bands) instead of the measured dark-matter halo 
masses. Finally, our method can also be extended to use 
redshift-space information to assess the effect of large-scale 
bulk motions and the small-scale finger-of-god effect on the 
resulting filaments. 

We will not attempt to find all the filamentary struc- 
tures in the simulation, only those filament segments gen- 
erated between haloes above a given mass threshold (node 
pairs). Therefore, smaller filaments associated to less mas- 
sive nodes will be missed if they are not in the path (or part) 
of the selected nodes. 



3.1 Density field 

In this paper we distinguish between two different defini- 
tions of density; (i) the standard dark matter density traced 
by the particles in the simulation which we call DM density; 
and (ii) the density given by the halo positions and their 
virial masses which we call the halo density. It is clear that 
the halo density contains little information about the mass 
and structure that lie beyond the virial radii of the haloes, 
but as we will show it is still an appropriate proxy of the DM 
density in the simulation. It is clear that halo positions and 
their masses (or in the observational case, galaxy positions 
and luminosities) allow a clear by-eye detection of filamen- 



tary features at large scales (jColberg. Krughoff fc Connolly I 
l2005h . 

In general, the density and density gradient field of a 
distribution of points can be o btained using VT, in a similar 
approach to that adopted by lAragon-Calvo et al. I l|2007bl ) 
where they compute th e density field using Delaunay Tes sel- 
lation Field Estimator l|Shaap fc van de Wevgaertl2000D . In 
the present study we make use of the neighbour information 
for all the haloes to trace the halo density field as well as 
to compute a fast proxy for the halo density gradient vector 
field. VT also allows us to obtain the immediate neighbours 
of each halo (or galaxy if applied to o bservational data). 
The Voronoi Tessellation l|Voronoi lll908h technique is one of 
the best adaptive methods to recover a precise density field 
from a discrete distribution of points, with clear advantage 
over the method used in Smoothed Partic le Hidrodynamic 
or ot her interpolation based techniques l|Pelupessv et al."l 
2003). We compute the VT for the halo distribution defining 
a cellular-like structure, where each halo is associated to a 
region (or voronoi cell) in which any point inside this region 
is nearest to that halo than to any other. 

This voronoi cell defines a volume which used along with 
the enclosed mass, defines a very precise and adaptive mea- 
sure of the density of the cell. In the case of point masses 
(such as when using the DM particle distribution), one can 
measure the exact enclosed mass in each voronoi cell, and 
therefore compute a very accurate DM density field. Instead, 
in this work we use the halo positions along with their mea- 
sured virial masses. The VT computation is done in the 
same way as for particles, but the halo virial mass does not 
account for all the enclosed mass in the voronoi cell, it only 
includes the mass out to the virial radius. For instance, in 
low density environments the halo-to-neighbour distance is 
much larger than the virial radius, and therefore the mass 
enclosed in the voronoi cell given by the virial mass of its 
central halo is underestimated. The opposite occurs in dense 
environments where the voronoi cell volume of a halo can be 
even smaller than their virial sphere due to close neighbours; 
in this case there is an overestimation of the enclosed mass 
in the voronoi cell. As this method does not require absolute 
density values but only the relative highest density path be- 
tween nodes (mainly given by the collapsed mass) the use 
of the halo density would increase the contrast of filaments 
improving the ability of the method to follow their high rel- 
ative density path to some degree. 

We argue that in the high density end, the halo density 
over-estimation is not important for our purposes since i) we 
will not consider subhaloes or haloes inside the virial radius 
of nodes (the most massive haloes) , ii) the inter halo distance 
becomes comparable to the virial radius at halo densities 
much greater than the average density along the filaments, 
and therefore only a few haloes considered in our analysis 
will suffer this overestimation. As a result, most of the haloes 
that will present an overestimated density will be nodes, 
and the remaining affected fraction will be located around 
nodes and in the central sections of the filaments, where 
their filament membership will be ensured, independently 
of the overestimation of their density. 

In low density regions the voronoi cells of haloes are al- 
ways much larger than their virial spheres which produces an 
underestimation of the density; later in this section we will 
work on diminishing this problem by using an approximation 
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assuming Navarro, Frenk & White (1997, NFW) profiles, to 
define the characteristic DM density between two haloes. 

Before moving on to the calculation of the character- 
istic density between haloes, we will analyse in more detail 
the differences between the halo and DM densities. For a 
smooth density field, such as is the case of fields traced by 
DM particles, the Hessian matrix can be computed with 
high accuracy to find the filament components easily. But 
the process is more complicated in the case of having only 
the positions of haloes and their virial masses. This is due 
to the sparse coverage of haloes, their variable masses, and 
the loss of information regarding the mass located beyond 
the virial radii of haloes. In order to understand the impor- 
tance of these issues we will look at the relation between 
average halo to neighbour separation (Dip) and its voronoi 
cell volume. 

In order to recover the real DM density field as best 
as possible using only halo positions, one needs to take into 
account that, 

• In high density environments the voronoi cell volume 
is related to the local mean inter-particle distance, i.e., the 
mean neighbour distance Dip. The left panel of Figure [5] 
shows a very tight relation between these two quantities for 
the full halo population. In the figure, the dashed line shows 



the V oc D\ p relation, whic h is very useful for halo detec- 
tion methods such as FOF (|Davis et al. |[l985l ). where the 
particle separation is used to connect particles above a given 
density threshold. In the case of having only halo positions, 
we find that this relation breaks down at lower densities (as 
can be seen in the left panel of Figure [2)| . The origin of this 
departure from the distance vs. volume relation is the com- 
plex shaped] developed by Voronoi cells at such low densities 
as result of greater standard deviation in the computation of 
Dip due to a low neighbour count and inter-halo distances 
falling within a wide range of possible values. Another pos- 
sibility is that shot noise is affecting our estimates, but this 
should not be the main source in our case since haloes mark 
the highest peaks in the density field, and we use a rel- 
atively large minimum number of particles per halo. This 
implies that the local dumpiness of a set of haloes in low 
density environments is only poorly related to its density; 
this may pose a challenge to the search for the backbone of 
filaments. This effect is negligible when obtaining the den- 
sity field using DM particles since these typically produce 

1 We refer as complex cell shapes to non-spherical or non- 
poly hedric like shapes, produced when having few neighbours at 
non-uniform distances. 
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Figure 2. Left panel: Voronoi cell volume vs. mean neighbour separation, Djp, for all the haloes in the simulation. Right panel: Voronoi 
halo Density vs. Djp for all the haloes. Dashed lines represent the relation V = Dj p (left panel); in the right panel it corresponds to 
V = median(My m) / p. 



a smoother spatial coverage and therefore a much smaller 
fraction of these will be surrounded by Voronoi cells with 
complex shapes. 

• The spread in the virial masses of the haloes, introduces 
a scatter in the relation between mean neighbour separation 
and the halo density (right panel of Figure(2]| with respect to 
the resulting relation from using only cell volumes. Therefore 
the halo density can only be used as a proxy for the matter 
density, and will serve to choose which halo pair will have the 
highest local DM density; by choosing the neighbour with 
the greatest halo density, it will probably be the nearest 
one and very likely the correct choice. This will sometimes 
not be true, for example when two or more neighbours have 
similar halo densities. Consider for example two neighbour 
cells with almost equal densities, but one having F times 
more mass and F times more volume than the other(F > 1); 
if we make a simple estimate of the DM density for the region 
lying between the halo and these two neighbours using NFW 
profiles for each halo, we will find that the path connecting to 
the smallest and closest neighbour will have the highest DM 
density. Later in this section we will apply this correction to 
our VT density estimates. 

We now estimate the local DM density between a halo 
and its neighbours, which we call the characteristic DM den- 
sity p*. As we have shown, the halo density estimate is rela- 
tive and it is only used to find the neighbour with the highest 
local DM density from all the possible halo-neighbour pairs. 
This density is an approximation that depends on the halo 
masses and inter-halo distances, and therefore it is probably 
safer not to compare it to the real DM density field given 
by the DM particles. Due to these considerations, in order 
to find the path of highest local DM density connecting two 
nodes, we need to add conditions on when and how to use 



of the halo density field. To estimate the characteristic DM 
density p* between the i-th halo and one of its neighbours, 
halo j, we will have two cases depending on the relation 
between their separation and their virial radii, 

i) Dij < Rvm{i) + Rvm{j) ■ p* = ki p(j), 

where Dij is the distance between haloes, p is the halo den- 
sity, and fci is a constant which includes the halo density 
of halo i common to all its neighbours. The fraction of 
halo pairs which satisfy this condition is very low and cor- 
respond to nodes and their immediate neighbours (haloes 
which are linked gravitationally); here the halo density is a 
good proxy for the DM density, and even a possible over- 
estimation of the halo density due to cell volumes smaller 
than virial spheres is positive for our purpose, since gravi- 
tationally linked haloes should have the first priority at the 
moment of choosing the halo-neighbour path to form the fil- 
ament skeleton. In this case, the segment connecting haloes 
i and j will have the maximum characteristic DM density 
among the other immediate neighbours. 

ii) > R V m(i) + Rvm(j) : p* = k 2 p(j) r)' 1 f(M h Mj). 

Most halo pairs fall in this second case. Here we use NFW 
profiles to estimate a proxy of the characteristic DM den- 
sity between two haloes. This proxy consists on the mini- 
mum DM density present in the path between two haloes, 
obtained by extending NFW profiles beyond the halo virial 
radii (this is a good approximation since the average of the 
inter-halo separation in the filament backbones is 4.80±0.03 
times the sum of the virial radii of the two neighbour haloes, 
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see Section 4.1). In the equation, Mi and Mj are the halo 
masses, the 77 factor represents the break-down of the rela- 
tion between inter-halo distance and voronoi cell volume, 
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and fci and are constants intended to provide the con- 
tinuity between both densities at r\ = 1, and Q = 1; 
M* = 10 12 ' 5 h _1 M Q is the constant in the iBullock et al. I 
concentration vs. mass relation. The rj parameter 
appears naturally in this approximation where its value is 
usually greater than one; therefore, two haloes with high 
masses and high voronoi halo densities will have lower p* if 
their separation is large, as can be the case in regions with 
a low number density of discrete points. 

The DM density between two haloes will be used as 
segment weights in the search for the path connecting two 
nodes, in a similar way to that used in the search for the 
shortest path in graph theory; therefore, the filament back- 
bone or skeleton is the result of solving for this graph, which 
has several different approaches in the literature (Biggs, 
Lloyd & Wilson, 1986). 



3.2 Input parameters 

We detect filaments using nodes above a fixed minimum 
mass. This choice is necessary since the filamentary structure 
is found at different scales; there are even fil aments inside 
filaments or inside clusters (|Bond et al. II1996I ). 

In addition to the minimum node mass, other param- 
eters will be necessary since otherwise it is always possible 
to find the highest density path connecting any two nodes. 
However, our aim is to involve only the lowest number of 
parameters possible, which include the following, 

• A minimum density threshold for the galaxies or haloes 
which form the backbone of a filament. This density refers to 
a minimum characteristic DM density (defined in the previ- 
ous subsection) along the consecutive halo pairs which form 
the filament backbone. There is no fixed physically moti- 
vated minimum value for this quantity, but we are interested 
in the filaments which are at least noticeably above the lo- 
cal background density, i.e. filament backbones above a few 
times the mean density. We will use this minimum density 
as a quality parameter for the detected filaments, since the 
higher this density for a filament is, the stronger the den- 
sity gradients and filament-like potential will be, with more 
haloes bounded to them. 

• A maximum gap threshold for the galaxies or haloes 
which define the backbone of the filament. A measure of 
the gaps in a filament is given by max(DsK / < Dsk >), 
the maximum distance divided by the average distance be- 
tween all pairs of consecutive skeleton members of an indi- 
vidual filament. Large values for this parameter imply large 
gaps between two filament sections. Gaps are an important 
problem, particularly for low density filaments. Again, this 
parameter will not define a limit on what is identified as 
a filament, but will be used as another quality parameter 
since the smaller this value is, the more continuous and uni- 



form the filament will be, with less noticeable gaps in the 
backbone. 

• After the definition of the backbone or skeleton of the 
filament has been completed, we select the members of the 
filament. This is done by analysing which neighbours are 
gravitationally linked to the filament and will collapse into 
the skeleton or remain within the filament for at least a 
given amount of time. We define a timescale £f, which is 
the maximum time allowed for the orbit of a halo in the 
plane perpendicular to the filament, assuming it is gravita- 
tionally bound (in this plane). Since the peculiar velocities 
of the haloes in the numerical simulation are known, we can 
calculate which haloes are bound to the filament; we use this 
information to characterise an average timescale and the as- 
sociated radius out to which bound haloes can be found. 
This will help to implement this filament identification in 
the case of observational data with no available information 
on peculiar velocities. 

It is complicated to define physically motivated density 
and gap thresholds for each filament analogous to the virial- 
isation density for the spherical collapse model. The reasons 
behind this are the complicated filament shapes and their 
continuous feeding of their node haloes or clusters. There- 
fore, we will use these parameters to assess the quality of a 
filament; filaments will be better defined if their minimum 
backbone densities are high and their largest gaps are small. 

The reasons behind the choice of these two parameters 
to define the quality of filaments are the following. A fila- 
ment is a region in the universe where the gravitational col- 
lapse of matter occurs mainly towards a line (continous but 
not necessarilly straight); therefore we have a cylindrical-like 
density profile with its associated cylindrical-like potential. 
Following this principle, and at the scales we are interested 
in in this paper (filaments between high mass haloes), we will 
assume a filament is of higher quality than another one if it 
is more likely to satisfy the previous conditions. A stronger 
cylindrical-like density profile (indicated by the DM den- 
sity between consecutive halo pairs in the skeleton) above 
the background will produce a stronger collapse of matter 
towards the skeleton, and smaller gaps between filament 
backbone members will better guarantee the continuity of 
the filament. The complex geometries and different scales 
characterising filaments, along with the facts that there is 
no known density profile a filament should follow and that 
they are unstable structures, make it difficult to set the val- 
ues for these two parameters that will ensure a high quality 
sample of filaments. Instead we simply assume that a higher 
characteristic density and smaller gap parameters imply a 
higher quality filament. 

3.3 Description of the algorithm 

Figure [3] shows a cartoon depiction of some of the steps fol- 
lowed by the algorithm to identify filaments for a particular 
node pair; in the figure, circles represent halo positions and 
their virial radii.We identify filaments in the following way, 

(i) We select a node tracer pair (indicated by blue circles 
in the figure). 

(ii) We follow the segments of highest local DM density 
given by the characteristic density p*. This defines the fila- 
ment backbone or skeleton. For this we define a set of thresh- 
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Figure 4. Four examples of detected filaments. The red solid lines show filament skeletons, the blue dashed lines show the re-centred 
skeleton. The white asterisks correspond to haloes at distances from the filament r < ro, whereas blue squares show haloes at distances 
r < r\. The red triangles show haloes with Ep < 0. 



old densities pth(i) with i — 1..N, in the range set by the 
minimum and maximum densities in the full density field. 

(iii) For each node we generate a list of neighbour haloes 
just outside the virial radius in the half hemisphere that 
points to the other node. These neighbours will be labeled 
as start haloes associated to the node from which we will 
start the filament search. End haloes will be the neighbours 
associated with the other node in the node pair in the half 
hemisphere pointing back to the start node. In panel a) of 
Figure O the blue dotted lines indicate the half hemispheres 
of the nodes that point to the other node; red circles mark 
the haloes at the start and end nodes. 

(iv) The first attempt at identifying a filament is done 
starting at the highest density threshold pthii = 1). 

(v) The process is iterative selecting the start halo with 
the highest local DM density with respect to the start node, 
characterised by a local density greater than pthip)- A halo 
that satisfies this condition becomes part of a possible skele- 
ton, and we search for neighbours of this new skeleton mem- 
ber using the same conditions. If there are no new neigh- 
bours satisfying this, we go back to the previous halo from 
where we will choose a different neighbour to restart the pro- 
cedure. Panel b) depicts this step. The colours of the lines 
(solid and dashed) connecting pairs of haloes correspond to 



the local characteristic DM density (densities are shown in 
the colour-scale bar at the bottom of the panels). As can 
be seen, we start with the maximum characteristic density 
threshold pthii — 1) denoted by a vertical black line in the 
colour bar. We choose the start halo (the one connecting 
with the start node located near the bottom of the panel) 
which has four neighbour candidates (connected by dashed 
lines to the start halo) for skeleton members, but two neigh- 
bouts are neglected since they are also start haloes. This 
leaves two remaining candidates, but none of them are char- 
acterised by densities higher than the threshold, and we are 
not able to find a filament at this density threshold. 

(vi) We repeat the last step with a different start halo 
until any of the end haloes of the other node are reached, or 
until there are no more haloes satisfying these rules. 

(vii) If no connection to the other node is found, we move 
down to the next lower density threshold step pthii + 1), 
and go back to step v. Panel c) shows the skeleton after 
lowering several times the density threshold down to the 
point where the skeleton contains four members (connected 
by the solid lines). However, the fourth skeleton member has 
no neighbour candidates (connected by dashed lines to the 
fourth member) with characteristic DM density greater than 
the current threshold. 



Automated detection of Filaments in the LSS 9 



(viii) We will always find a set of connected points (a fil- 
ament backbone) between two nodes for a sufficiently low 
value of pth density. Higher values of this density imply 
stronger filament contrasts. Panel d), shows the result when 
a first skeleton was completed between the two nodes, for a 
sufficiently low density threshold. 

(ix) We re-centre the local centre of mass of the filament 
skeleton using its immediate Voronoi neighbours. 

Having a well defined backbone, we start adding skele- 
ton neighbours to the filament and computing filament char- 
acteristics, in the following way, 

(i) For any given halo k we find the nearest skeleton mem- 
ber j (shown in panel e) 

(ii) We measure the mass contained in a cylinder around 
the skeleton at the position of the skeleton halo j. The 
cylinder height is H = [Dj^+j + Dj t j-i)/2 and its radius 
R — Dk,j .Using this mass and the difference between the av- 
erage velocities of the haloes within that cylinder and that 
of halo k, projected in the plane perpendicular to the cylin- 
der, we compute the total halo energy in the plane, Ep. In 
panel e) of Figure [3] the cylinder is depicted by black dashed 
lines. The cylinder axis(middle black dashed line) is tangent 
to the filament at halo j as inferred using the two immediate 
neighbour skeleton members. 

(iii) We compute the orbit time t around the cylinder for 
halo k assuming that the distance D^j is the semi-major 
axis of the orbit. This timescale only uses information on 
the potential energy and does not require peculiar velocity 
data. 

(iv) We select all haloes with Ep < and calculate their 
median orbit time t\\ we define n as the radius containing 
80% of these haloes. This sample can only be obtained from 
haloes with peculiar velocity information. 

(v) We select all haloes with Ep < and t ^ tj?, with 
ip a fixed input parameter, and we define ro as the radius 
where 80% of these haloes are contained. This defines a sam- 
ple using Ep measurements and it therefore needs peculiar 
velocity information to be constructed. 

(vi) We select all haloes with t ^ tp, and define r-2 as 
the radius where 80% of these haloes are contained. This 
selection can be done with position and mass information 
alone and does not require dynamical information. 

(vii) Finally, we also select all haloes with t ^ ti, and 
we define r3 as the radius containing 80% of these haloes. 
This selection also requires velocity information and is used 
to assess the importance of the binding energy condition 
against that of the orbital timescales. 

All haloes closer to the skeleton than ri will be se- 
lected as filament members in the simulation. Panel f) of 
Figure [3] shows the resulting filament, where blue circles 
correspond to haloes belonging to the new filament; the re- 
maining nearby haloes are too far away from the filament 
and do not satisfy the membership conditions. 



4 RESULTS 

Figure [4] shows four detected filaments in the simulation, 
where the halo density projected onto the x — y plane is 
shown in a colour scale, the skeleton is shown as red lines, 
and the re-centred skeleton as blue dashed lines. The nodes 



are indicated by circles with radii equal to the halo virial 
radius. White points denote all haloes lying closer than ro 
from the filament skeleton, and blue boxes denote haloes 
closer than n. The red triangles are for haloes with Ep < 0. 
All the filaments contain segments with only either a few or 
no bound haloes, at least according to our definition. 

We bear in mind the possibility of undetected bound 
haloes since in our energy calculation we do not take into 
account nearby structures other than the filament. In order 
to produce a more precise energy calculation one would need 
to use velocities from other sections of the skeleton instead 
of only from the nearest skeleton section; filaments show 
a very complex velocity structure where nodes sometimes 
move towards each other (they may merge in the future) 
or away from each other, making filaments suffer stretch- 
ing, elongations, torsions, and even rotations. However, the 
incompleteness in the sample of bound haloes should not af- 
fect our estimate of the mean effective radius of the filament 
(ri) which we use to define filament membership. 

In the upper-left and bottom- left panels of Figure [4] the 
filaments show excellent density contrasts, but also show a 
gap (near the top node in the upper-left panel, and near 
the left node in the bottom- left panel). This shows the im- 
portance of adopting a gap parameter that allows the exis- 
tence of these features in selected filaments to some degree. 
The filaments in the right panels are of higher quality than 
those on the left since they do not show important gaps. The 
section of the filament on the upper-right panel seems not 
to follow the highest density path due to projection effects 
(the filament follows a path that enters the page, along the 
z-axis) . 

4.1 Filament properties 

We apply the method to the numerical simulation described 
in Section 2, using a minimum skeleton characteristic density 
p*min = Spmean and no gap restriction, limiting the node 
pairs to relative distances lower than 65Mpc/h. 

Out of the 3385 node pairs, 1326 are successfully con- 
nected via filaments; we will refer to this first identifica- 
tion as the full sample. We select an additional subsam- 
ple of 467 filaments which satisfy the additional condi- 
tions of p*min above the median of the full sample, and 
max(DsK)/ < Dsk > below the median; this sample is 
termed the high-quality subsample and contains 33% of the 
filaments in the full sample. The separation between back- 
bone members in the full and high-quality samples are, on 
average, 4.80 ± 0.03 and 4.19 ± 0.03 times the sum of their 
virial radii, respectively. As was mentioned above, all the 
detected filaments connect nodes separated by at least the 
sum of their virial radii. Figure [S]shows the relation between 
gap and density parameters for the detected filaments which 
show clear trends of larger gaps at lower densities. 

Figure [5] shows the dependence of the quality param- 
eters on node separation for the full sample. There are 
clear correlations, particularly for filaments shorter than 
20Mpc/h, which suggests that shorter node separations pro- 
duce higher quality filaments. 

When studying the properties of the filaments detected 
using our automated procedure, it will be useful to com- 
pare with a pr evious detection. In particular, we wi ll use 
the results from lColberg. Krughoff fc Connolly I (|2005T 1 who 
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Figure 6. Gap (left) and density (right) quality parameters as a function of node separation. 
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Figure 5. Filament quality parameters. Minimum skeleton den- 
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detected 228 filament in a DM simulation by eye using the 
smoothed DM density distribution. This filament sample can 
not be compared directly with our results, since the selection 
criteria are very different. However, both samples are the re- 
sult of restricting the search to filaments connecting neigh- 
bouring haloes above 10 1i _1 Mq. The main differences be- 
tween the two samples arise from, i) Colberg et al. use the 
distribution of DM particles whereas we use halo positions, 
ii) they look for filaments using the 12 nearest haloes inside 
cylinders of 7.5h~ Mpc of radius aligned along the node- 
node axis; in our case we look at all possible neighbour node 
pairs given by the voronoi tessellation with no volume con- 
strain, iii) Colberg et al. define a true detection based on a 



visual criterion instead of using quality parameters, iv) they 
discard node pair connections when other clusters lie inside 
the innermost 5h _1 Mpc from the node-node axis, and we 
discard node pair connections when another cluster is closer 
than 2 times its virial radius to the filament skeleton, v) 
they divide their sample in straight, off-centre and warped 
filaments. Therefore, the reader must bear in mind that com- 
parisons between these two samples, are not intended to val- 
idate any of the two samples, but to find general filament 
properties which are less sensitive to different selection cri- 
teria. 

The node pair connections are given by the voronoi tes- 
sellation method, which instead of selecting the n nearest 
neighbours, chooses neighbours such that the line that con- 
nects the pair passes only through the voronoi cell around 
each node. This ensures that any point along the segment 
is nearest to one of the two nodes and not to other haloes. 
The node pair count of the full sample as function of the 
node separation is shown in Figure [7J as an orange dashed 
line (the scale of the counts in 5h _1 Mpc bins is given by 
the right y-axis). The number of pairs grows almost linearly 
with the separation almost up to 40h _1 Mpc, and then it 
decreases for larger node separations. In addition, Figure [7] 
shows the fraction of node pairs with detected filaments as 
a function of node separation (left y-axis scale). The full 
sample (solid black bars) is characterised by a decreasing 
fraction of connected pairs via filaments as the separation 
increases; this fraction is nearly 90% for separations shorter 
than 5h _1 Mpc, and at the largest separations the fraction 
is reduced to 30%. In the case of the high quality subsample 
(red bars) the abundance of filaments decreases much faster 
with fractions below 25% for nodes separated by more than 
20h _1 Mpc. 

Figure [7] also shows the fractio n al abu ndance obtained 
by IColbergTKrughoff fc Connolly I |2005l ). Even though 
their selection procedure is different from ours, the resulting 
dependence of this fraction with pair separation is similar to 
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Figure 7. Fraction of node pairs with detected filaments as 
function of node pair separation (left y-axis). The green bar on 
the first bin shows the fraction of node pair connections separated 
by less than the sum of their virial radii. The total node pair count 
as a function of node separation for the full sample is shown as 
an orange dashed line (its scale is indicated on the right y-axis). 
Black bars show the same fraction for the full sample of quasars, 
red bars are for the high quality subsample, and the blue bars are 
for the Colberg et al. sample detected by eye. 



our results for the high quality subsample. Our method does 
not consider haloes inside the virial radii of nodes, which 
means that we do not detect most of the filamentary struc- 
ture connecting two nodes separated by distances shorter 
than the sum of their virial radii. The green bar shown for 
separations shorter than 5h _1 Mpc in Figure [7] indicates the 
fraction of node pairs whose separation is shorter than the 
sum of their virial radii for this range of separations. Most of 
the pairs rep resented by the green bar should be connected 
by filaments (|Pimbblet et al. 112004 ) . since th ese are overlap- 
ping bound systems which share matter (i.e. iDietrich et ahl 
120051 : iTittlev fc Henriksen Il200ll ): furthermore, this behav- 
ior should extend up to node separations of a few virial radii 
(approximately three times the virial ra dius), which can be 
associated to the infall region of halo es (|Diaferio fc Gellerl 
ll997l ; |Pivato. Padilla fc Lambas1l2006l ). Taking into account 
the mass resolution of our numerical simulation, the mass in 
such bridges is mostly in the form of a smooth DM particle 
distribution, with only a few subhaloes aligned within the 
bridge. This makes it more difficult to detect them with our 
method even if we also used subhalo positions; as a conse- 
quence we have chosen not to include them in the search. A 
possible way to overcome this would be to use the DM par- 
ticle distribution, or to run re-simulations of these regions 
with higher resolution, enough to resolve several subhaloes 
per node. The result of such a study would likely change 
our fraction of detected filaments for node separations be- 
low 5h -1 Mpc, which we are underestimating at present; in 
the case o f lColberg, KrughofT fc Connolly I (|2005i ), they find 



that most of the halo pairs within this range of separations 
are connected via filaments. 

In the case of fi laments detected in the 2dFGRS, 
iPimbblet et~ai~l (|2004l ) find a fractional abundance of fila- 
ments similar to our full sample results; however, their selec- 
tion criteria are also different from the one we have applied 
to the simulation. In particular, they also identify filaments 
by eye and use galaxy positions; therefore in order to make 
an appropriate comparison it would be necessary to apply 
our method to realistic 2dFGRS mock catalogues, or directly 
on the 2dFGRS catalogue. 

The main properties of the detected filaments are shown 
in Figure [5] The top- left panel shows the distributions of ti 
(the median orbit time for haloes with Ep < 0), where it 
can be seen that the high quality filaments are characterised 
by lower orbit times as expected since these filaments have 
higher density contrasts and are more concentrated than the 
full sample. The samples shown in the figure are obtained by 
setting tF = 2to (vertical red dashed line) which is slightly 
lower than the median of t\ (indicated by the vertical blue 
dashed line). This latter value can be used when detect- 
ing filaments without dynamical data since the orbit time 
distributions shown here are relatively narrow (most of the 
filaments show similar orbital timescales). 

The top-right panel of the figure shows the distribu- 
tions of the parameters ro and ri (line types are indicated 
in the figure key) described in the previous section. As can 
be seen, a fixed orbit time produces a narrow distribution 
of ro but a wider distribution of r\ which is obtained using 
t\. However, the peaks of both distributions are located at 
approximately 1.3Mpc/h. It is also noticeable a very slight 
shift towards smaller radii for the high quality subsample in 
both cases, an effect which is stronger for the r\ parame- 
ter, indicating a dependence of the tj value with the quality 
of the filaments. Therefore, better quality filaments seem to 
be more concentrated while preserving similar thicknesses 
with respect to the filaments in the full sample. In addi- 
tio n, the figure also shows the scale ra dius r s computed 
by IColberg. Krughoff fc Connolly I l|2005l ) for their sample 
of filaments (blue bars). In their notation r s defines the ra- 
dius where the density profiles of straight filaments starts to 
follow a r~ 2 relation. Our definition of ri indicates a scale 
radius containing 80% of the bound haloes with orbit times 
below the median. Even though both definitions are con- 
ceptually different, they account for the scale radius where 
f» 50 — 80% of the filament mass is contained. In general, 
for a given filament, r s is a more precise computation of the 
edge of the filament, but requires the DM particle distribu- 
tion to be calculated; r\ is easier to compute since it only 
requires halo positions; however, it can underestimate the 
filament edges depending on the density profile and density 
contrast. Therefore, despite the fact that the comparison 
is made among two quantities with different definitions, as 
well as different filament samples, it is interesting to note 
that the distributions of r\ and r s show similarities; the lat- 
ter only shows a slight shift towards larger radii. As can 
be seen, the characteristic radius which defines a filament 
shows a narrow distribution with preferred values of 1 to 
2h _1 Mpc, even when using filaments of different quality or 
using a sample of filaments selected by eye. In all cases, how- 
ever, the lengths of the filaments are similar and are traced 
by halo nodes with masses above 10 14 h _1 M Q . 
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Figure 8. Distribution functions of properties of the detected filaments. The sa mples of filaments detected using our au tomated method 
are shown in different line types (explained in the key). The statistics from the IColberg. Krughoff fe Connolly I {2005) filament sample 
are shown as barred histograms. 



The Bottom-left panel of Figure [S] shows the distribu- 
tion of mass for different filament components (line types are 
show in the figure key); all the distributions are shown for 
the full sample of filaments. As can be seen, this tracer node 
mass selection produces skeletons and filament envelopes less 
massive than the filament nodes. Both of these two compo- 
nents show similar distributions, with differences only at the 
low-mass-end. Notice that when using either ro or ri the re- 
sulting filament mass is practically the same. This shows 
that the detection of filaments using a fixed orbit time ( 
when no dynamical information is available) will provide re- 
liable filament mass measurements. In the case of the high 
quality filaments, we find that the masses of the skeleton 
and the surrounding filament shells are lower than for the 
complete filament sample, since the former are shorter in 
length (as can be seen in the bottom-right panel of the fig- 



ure). We find no clear dependence of filament mass on their 
node masses. 

The bottom-right panel of the figure shows the distri- 
butions of node pair separation and of filament extension 
(line types are indicated in the figure). The filament exten- 
sion is obtained by adding the distances between consec- 
utive filament member positions (i.e. in a discrete number 
of segments) along the filament. The node separation is on 
average smaller than the filament length, which indicates 
that most filaments are warped. The distribution of node 
pair separation peaks at « 32Mpc/h for the full sample, 
and at m 15Mpc/h for the high quality subsample. The fil- 
ament lengths also show a peak at shorter values for the 
high quality subsample. When analysing the ratio between 
these two quantities in both, the full and high quality sam- 
ples, it can be seen that regardless of quality, longer fila- 
ments are more warped than shorter filaments; i.e. in the full 
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sample, filaments with node separations below 30h~ 1 Mpc 
are on average 13% larger than their node separation; this 
value increases to 40% for larger node separations. The 
filaments studied with Shapefinders in the Las Cam panas 
Redshift Survey (|Bharadwai. Bhavsar fc Shethl 12004! ) are 
characterised by lengths of 50h~ 1 Mpc to 80h _1 Mpc. In this 
work, we have found shorter high quality filaments but we 
have also required node pair separations < 65h~ Mpc. It 
should also be borne in mind that in most cases these fil- 
aments are only segments of considerably longer structures 
with more than two nodes (shapefinders are insensitive to 
the number of nodes in a filament). 

Figure [9] shows the relation between filament thick- 
ness (ri) and filament length. The error bars correspond 
to the standard deviation in the measurement of the me- 
dian of n, and are computed using the jackknife method. 
In the high quality subsample, we do not include fila- 
ments longer than 80h _1 Mpc due to low filament counts 
(< 10). As can be seen, there is a trend of thicker fil- 
aments for longer filament lengths in both samples (full 
and high quality). For the high quality sample, the me- 
dian value of r\ for filaments with lengths between and 
10h _1 Mpc is 1.11 ± 0.19h _1 Mpc, and for lengths between 
60 and TOh^Mpc it is 2.01 ± 0.29h _1 Mpc (a significance 
of more than 3rr for a difference between the longest and 
shortest filament lengths). This dependence can be a con- 
sequence of any or several of the following effects, i) all 
filaments feed their node haloes and shorter, less massive 
filaments will exhaust their mass first due to the higher in- 
fall velocity and node halo influence over a larger percent- 
age of the filament length (the influence ca n extend out to 
several virial radii. iDiaferio fc Geller Ill997l ). ii) shorter fila- 
ments are straighter than longer ones; therefore, in longer, 
warped filaments concave zones along the skeleton could at- 
tract haloes from larger distances, an effect that would be 
absent in straight-line filaments. The detailed study of this 
possibility is beyond the scope of this paper and will be 
treated in a forthcoming paper on filament shapes and en- 
vironments, iii) A higher probability to spuriously assign 
bound haloes at larger distances from the skeleton for longer 
filaments, but this is less likely since this effect is also present 
when using ro (which does not depend on a computation of 
energy) as a thickness indicator. 

We study the variation of the mass density along the fil- 
ament skeletons. Figure[TU]shows the average over-density as 
a function of the normalised node pair separation. It should 
be borne in mind that as we use the interpolated voronoi 
density obtained from the halo positions and their viral 
masses, the density only includes a fraction of the total mat- 
ter (DM particles beyond the virial radii of haloes are not in- 
cluded in this estimate). We exclude filaments with skeletons 
containing less than 6 haloes, and the figure only shows half 
of the filament length since the profiles are symmetrical (on 
average) . The figure shows a similar dens i ty pro file to those 
found bv lColberg. Krughoff fc Connolly I l|2005l ). where the 
over-density rises towards node centres, indicating that on 
average the infall regions of filaments extent up to 20% of 
the filament length. At larger distances from the nodes, the 
overdensity remains at nearly constant values of a few times 
the average density. The high quality subsample shows a 
similar profile although with higher density contrasts than 
the full sample. 
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Figure 9. Filament thickness (as measured by ri) as a function 
of filament length, for the full sample (black) and high quality 
subsample (red). 
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Figure 10. Average longitudinal filament over-density profile 
obtained using the interpolated voronoi density along the skele- 
ton, as function of the normalised node pair separation. We only 
show half of the filament length since the profiles are symmetrical, 
on average. 



We now study the number of filaments connected 
to individual nodes, and how this depends on the node 
properties. Figure [TT] shows the fraction of filaments con- 
nected to 0, 1, 2, ... filaments for the full sample (black solid 
lines), the high quality subsample (red solid lines), and the 
IColberg. Krughoff fc Connolly f (|2005h results (blue bars). 
The Poisson error amplitudes are shown as dashed lines for 
the full and high quality samples. In the full sample, most 
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Figure 11. Fraction of nodes connected to N filaments for the 
full sample of filaments (solid black lines), and for the high qual- 
ity subsample (solid red lines). In both cases, Poisson errors are 
shown by the dashed lines. The results from Colberg et al. (2005), 
are indicated as blue bars. 



nodes are connected to 4 — 6 filaments, indicating that allow- 
ing in all the detected filaments without applying any qual- 
ity constraints does not provide realistic results, bearing in 
mind the observational (Pimblett et al., 2005) and numerical 
simulation (Colberg et al., 2008) results on this statistics. A 
better agreement with these estimates is obtained when us- 
ing the high-quality subsample, in which case most nodes are 
connected to 2 filame nts (and the distribution is very sim- 
ilar to that from the IColberg. Krughoff fc Connolly 1 120051 
filaments). In general the number of filaments per node is 
strongly dependent on the quality of the filaments consid- 
ered; similar quality thresholds are needed in order to make 
meaningful comparisons. Given that the number density of 
filaments is three times higher for the full sample than for 
the high-quality sample (simply due to the total number of 
objects in each sample) it can be expected that the distri- 
bution of filament connections per node will also be a factor 
of three higher for the full sample, that is ~ 6 compared 
to ~ 2 connections for the full and high-quality samples, 
respectively (as ~ 83% of the nodes of the full sample of 
filaments are connected by high-quality filaments). 

Figure [12] shows the average number of filaments per 
node as a function of node mass. In all cases this number 
increases with the node mass. Errors, shown as dashed lines 
for the full and high quality samples, are obtained using 
the jackknife method; errors are not shown for the highest 
mass bin M > 1O 15 M0 (cyan hatched region) due to the low 
number of nodes (10) at this end. Nodes in the high qual- 
ity subsample are connected to an average of 1.87 ± 0.18 
filaments for the lowest mass bin, a number that increases 
to 2.49 ± 0.28 for M ~ 1O 14 ' 9 M ; the significance of this 
trend is higher than a 3er l evel. This behav i or wa s also 
observed in the 2dFGRS by iPimbblet et al. I (I2004TI, and 
in a numerical simulation (|Colberg. Krughoff fc Connolly I 



2005), clearly indicating that more massive haloes are more 
likely to have a larger number of connected filaments. This 
can be associated to the higher amplitude of clustering of 
more massive haloes characterising random gaussian fluct u- 
ation fields in a AC DM cosmology jPimbblet et al. II2004T I. 

There is a number of possible issues that could affect 
this statistics that need to be borne in mind, i) we do not 
use subhaloes, and therefore node pairs closer than the sum 
of their virial radii could present filaments which we do not 
detect. Such close pairs will be more abundant for more mas- 
sive haloes due to their higher local overdensities, therefore 
these undetected filaments could populate the high mass end 
of the figure [121 ii) To avoid repeated filament segments, we 
disc ard filaments which are closer than 2r v j r to a third node, 
and lColberg. Krughoff fc Connolly I (|2005l ) use a fixed value 
of 5h _1 Mpc for a similar proximity condition. In both cases 
we could be missing short filaments in dense environments 
where nodes are more massive, have larger virial radii and 
are more strongly clustered; in such places this proximity 
constrain could be excessive. In order to test this issue, we 
make a subsample of filaments applying the quality con- 
straints used for the high quality subsample, but allowing 
filaments closer to a third node when, a) the node pair sepa- 
ration is less than 10h _1 Mpc, b) the minimum density along 
the filament is greater than 10 times the mean density, c) the 
sum of the virial radii of the nodes is > 2.5h _1 Mpc, d) the 
filaments are close to straight-line shapes. These modifica- 
tions, in conjunction with the intrinsic properties of voronoi 
tessellations for the node pair selection, ensures that it is 
very unlikely that the short filaments in this new sample 
are repeated segments of other detected filaments. The rea- 
son behind this is that for larger node pair separations, there 
will be larger distances from a node to node axis to a third 
node. Otherwise the constrain of a common facet between 
node pair voronoi cells would not be fulfilled. This test sub- 
sample is shown as green long dashed line in Figure Q2] as 
can be noticed the relation of filament connections as a func- 
tion of mass becomes stronger. 

4.2 Application to observational data 

In the case of applying this method to galaxies, we can use 
luminosities instead of halo masses and detect filaments fol- 
lowing the path of highest luminosity density. In this 
the light of a galaxy is more concentrated than the mass it 
is safer to assume that the voronoi density traces that of the 
luminosity in both, the high and low density regimes. In this 
case, the filament quality can be defined using a luminosity 
density parameter as well as a gap parameter. However, it 
would become more difficult to measure a filament thickness 
since in general there would be no information on the mass 
and, additionally, there seldom is dynamical information to 
calculate binding energy conditions in galaxy samples. 

A possible way to apply this method could use the skele- 
ton brightness and a brightness threshold for filament mem- 
bership where i) the distribution of filament thickness, and 
ii) the relation between filament thickness and lenght match 
the results from a DM simulation where the settings on the 
quality parameters result in similar number densities of fila- 
ments. These tests, and an application to observational data 
from the SDSS are part of a forthcoming paper. 

In the case of a sample with estimates of galaxy masses 
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Figure 12. Number of filament connections per node as func- 
tion of node mass. Different line types correspond to sam- 
ples selected in this paper (identified in the figure key); the 
barred histogram corresponds to th e sample of filaments in 
IColberg. Krughoff fc Connolly I i2005T ). The hatched area shows 
the range of masses containing only 10 node pairs in our numerical 
simulation. 



but no dynamical data, such as in nearby galaxies, it would 
be possible to select filament members assuming that galax- 
ies are bound to the filament, and requiring orbit times lower 
than tF- In the simulation, as can be seen in Figure |8j using 
a fixed orbit time allows to recover a distribution of ro (see 
Section 3.3 for the definitions of ro, ri, ri and r%) which, al- 
though slightly narrower, peaks at the same radius as when 
using the full energy calculation. Also, the recovery of the 
filament mass is only mildly affected by the use of ro or ri 
to select filament members. 

Figure [Tol shows the relation between ro and r\. As can 
be seen, there is a linear relation between these quantities for 
n < median(ro) . Filaments in the high quality subsample 
show a very similar median ro and a slightly lower median n 
than the full sample, an effect that probably arises from the 
fact that filaments in the high quality subsample are shorter 
than in the full sample (see Fig. [HJ. In the case of the obser- 
vational data with masses, but no dynamical information, 
the method would only provide measurements of r2 which, 
when comparing the vertical long dashed and dotted lines 
in both panels, can be seen to provide a good approxima- 
tion to ro. As the relation between ro and n is reliable for 
thin filaments, r2 < 1.2Mpc/h, thick filaments will proba- 
bly suffer from an under-estimation of their real thickness, 
particularly if their quality is low. Regarding rz (horizontal 
dotted lines), it can be seen that their median values are 
very similar to that of ri, indicating that if one can esti- 
mate the collapse time of bound objects to the filament, the 
membership obtained using this estimated time will provide 
a good membership criterion. 



5 CONCLUSIONS 

We presented an automated method to detect filaments in 
cosmological simulations, using haloes above a fixed mass 
as tracers of filament nodes. In addition, we proposed possi- 
ble directions to improve this method to allow its use with 
observational data. As filaments cannot be treated as viri- 
alised structures as in the case of haloes, and as they are 
characterised by a wide range of lengths, it is a difficult 
task to identify them automatically. As a result these have 
been mostly identified by eye. In this work we detect fila- 
ments using an automated algorithm that provides two fila- 
ment quality parameters, i) a minimum skeleton character- 
istic density, and ii) a gap parameter given by the maximum 
distance between consecutive skeleton neighbours divided by 
the average consecutive skeleton neighbour distance in indi- 
vidual filaments. A small gap parameter and a high density 
parameter, ensure the best quality for a filament. The latter 
condition is equivalent to request a high density contrast. 

In our method we define the width of filaments using 
the median radius (ri ) that contains the haloes gravitation- 
ally bound to the filament in the plane perpendicular to 
the filament skeleton, and that are characterised by orbit 
or collapse times below an upper threshold. An application 
of the method to data without dynamical information can 
be done since the radius r\ shows a good correlation with 
ro and r2 (ri is obtained assuming that all the galaxies are 
bound to the filament and computing their orbit times based 
only on their positions and masses); the members are then 
selected requiring orbit times below a fixed time tF- The re- 
lation between ri and r2 is one-to-one for thin filaments be- 
low ro ~ 1.2Mpc/h; in thicker filaments r2 tends to slightly 
under-estimate the actual width of a filament. 

We have presented several filament properties which can 
be studied in observational catalogues such as the SDSS. In 
particular, a subsample comprising the 33% highest quality 
filaments in our numerical simulations shows very similar 
properti es to filaments detected by eye in num erical simula- 
tions bv lColberg. Krughoff fc Connolly I (|2005l ), 

• Filament lengths are mostly concentrated below 
50h _1 Mpc, but can extend to up to 150h _1 Mpc 

• Shorter filaments are characterised by more straight- 
line geometries than longer filaments. Filaments with node 
separations below 30h _1 Mpc are 13% longer than the dis- 
tance between their nodes; this increases to 40% for larger 
node separations. 

• The distribution of filament widths is relatively narrow 
and shows a clear peak at d = 3h _1 Mpc. There are indica- 
tions of an increase in the filament thickness as the filament 
length increases. 

• Nodes are connected on average to 2 filaments, this 
number increases slightly with the node mass, reaching ~ 3 
filaments per node for masses close to 10 15 Afo 

• In the infall region around nodes the average central 
skeleton density can be as high as a hundred times the mean 
density; at larger distances the density drops to a few times 
the mean density, and maintains a roughly constant value 
along 20 — 80% of the filament length. 

• There is a strong relation between length, quality, and 
straightness in the filament shape, where shorter filaments 
have better quality and are closer to straight-line geometries. 
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Figure 13. ro vs. ri for the full sample of filaments (left panel) and for the high quality subsample (right). The vertical and horizontal 
dotted lines show the median values of ro and r\ (respectively), which are quantities obtained using the full binding energy calculation. 
The dashed lines show the median values of T2 and r% (vertical and horizontal lines, respectively), which are the equivalent to ro and r\ 
for the case with no dynamical information (and therefore no energy calculation). 



Similarities of the high-quality sample with the 



IColberg. Krughoff fc Connolly I (|2005l ) results seem to in- 
dicate that the natural by-eye criteria are strongly re- 
lated to our quality parameters; a detection by eye se- 
lects high density contrasts and few gaps. We stress the 
fact that did not intend to match the properties of the 
IColberg. Krughoff fc Connolly I (|2005f ) filaments, instead we 
simply chose the mean values of minimum density and gap 
parameters to define our high-quality sample. 

The filament properties we have studied in this work are 
focused on the general characteristics of filaments. There 
remain many specific properties of filaments and of their 
galaxy populations which can be related to several recent 
results such as (i) the halo clusterin g dependence on the halo 
mass and on its formation time l|Gao. Springel fc White I 
2005), (ii) the correlations betwe en halo concentration an d 
spin with the local environment l|Avila- Reese et al. II2005I ). 
(iii) the fa ct that ga l axy s pins are strongly aligned along 
filaments l|Pimbblet I l200ot ). ( iv) the results using semi- 
analytic models obtained by iGonzalez fc Padilla I ([2009) 
which show several variations of galaxy properties with 
the local and large scale environment, as well as (v) other 
results showing that galaxy formation should be strongly 
dependent on the large scale environment starting from 
their early stages of development, due for example to the 
delayed reionisation of filaments with respect to clusters 
as shown by hydro-si mulations of the intracluster medium 
ijFinlator et al. II2009T ). A first step will be to compare ob- 
servational galaxy properties in filaments, in particular their 
colours, star-formation rates and luminosities with results 
from semi-analytic models, to characterise some of the pre- 
viously mentioned environment effects. 

Several studies of galaxy properties in clusters and voids 
have opened the possibility to expect important variations 
in the properties of haloes or galaxies while embedded in 



filament-like environments, since the populations of galaxies 
and haloes are very different in voids and clusters. By con- 
verging to a standard filament classification and detection 
method, the study of galaxy properties and halo assembly 
in filaments can be carried out with great detail to help 
understand the reasons behind these important population 
changes. 
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